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New method is developed for calculation of single-particle 
occupation numbers in finite Fermi systems of interacting par- 
ticles. It is more accurate than the canonical distribution 
method and gives the Fermi-Dirac distribution in the limit of 
large number of particles. It is shown that statistical effects 
of the interaction are absorbed by an increase of the effective 
temperature. Criteria for quantum chaos and statistical equi- 
librium are considered. All results are confirmed by numerical 
experiments in the two-body random interaction model. 

PACS numbers: 05.45. +b, 31.25.-v, 31.50.+W, 32.30.-r 



Studying numerous problems related to many-body com- 
pound states one needs to know whether the laws of statis- 
tical physics can be used in the description of a particular 
quantum system with finite number of interacting particles 
(nucleus, atom, etc). The aim of the present work is to de- 
velop an accurate method for calculation of occupation num- 
bers (and other statistical average values) in finite isolated 
systems and compare the results with the Fermi-Dirac (FD) 
distribution. The latter is a priori valid in infinite systems 
of non-interacting particles. We show that the accuracy of 
the FD distribution in finite systems of interacting particles 
can be improved by introduction of the effective temperature 
which absorbs statistical effects of the interaction. 

For finite systems of interacting Fermi-particles the occu- 
pation numbers n 3 of single-particle orbitals can be found if 
we know the expansion of exact eigenstates \i) in terms of 
Slater determinants ( "shell model states" ) 



*>=£k 



(k\h s \k) 



i) = Y,C^\k); \k) =a\...ai |0) 



(1) 



(2) 



Here, h B = a\a a is the occupation number operator. 

If the residual interaction between the particles is strong 
enough, the expansion (j^) can be treated as a "chaotic" 
superposition of the basis states \k) . In this case mean 
square values of the expansion coefficients are smooth func- 
tions of the energy difference between the energy Ek of the 
basis state \k) and the energy E^' of the exact state 
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In numerical studies of the Ce atom |fj, the s — d nuclear 
shell model Q and random two-body interaction model 
ji| it was demonstrated that typical shape of exact eigenstates 
( "spreading function" ) F^ practically does not depend on a 
particular system and has a universal form characterized by 



the spreading width T . For example, for Ce atom Jl] the 
squared Breit-Wigner shape of F^ 1 ' has been found in a good 
agreement with numerical data, F^ tx [(E^ — E^ — A^') 2 + 
r 2 /4]~ 2 , where A' 1 ' < T is a small shift (see below) which 
in the zero approximation can be neglected. The value of 
r 1 1/ 2 

r = 2 
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can be found using the following relations 



for the spreading width of basis components, 
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with Hk P standing for non-diagonal Hamiltonian matrix ele- 
ments defined by residual interaction V. For example, in the 
model of n particles distributed over m orbitals we have || 
(r/2) 2 = (AE) 2 = Vgn(n - l)(m - n)(3 + m - ra)/4. Here 

V 2 = \V at | is the mean square value of two-body residual 
interaction matrix elements. 

Thus, we can rewrite Eq. (|l|) in the form of the " 
F— distribution" which gives the actual distribution of occu- 
pation numbers in finite Fermi systems, 
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where ni fc ' = (k\n B \k) equal or 1, e s are the energies of 
single-particle orbitals, and u sp = V S p^ S p is the diagonal ma- 
trix element of residual interaction. In practice, the second 
term in the definition of Ek can be substantially reduced by 
an appropriate choice of the mean field. The denominator in 
([|) stands for the normalization of the F— distribution. This 
"microcanonical" distribution is convinient for numerical cal- 
culations. 

It is instructive to compare the F— distribution (^) with 
occupation numbers obtained by making use of the standard 
canonical distribution, 
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were T is the temperature. The difference between Eqs. (^) 
and (JsJ) is that the summation in Eq. (^J) is performed over 
simple basis states while in the canonical distribution the 
summation is carried out over exact eigenstates. Another 
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difference is that in Eq. (^) the occupation numbers are cal- 
culated for a specific energy E of system unlike specific tem- 
perature T in Eq. (|5[) . However, Eqs. and (^) can be 
compared with each other by setting E — {E) T . 

It is important to demonstrate that the F— distribution (^) 
tends to the standard FD distribution in the limit of large 
number of particles. By performing explicitly the summation 
over n a = and n a — 1, the expression (Q) can be written in 
the form 



n 3 (E) 



+ Z a (n- l,E-t a ) 
Z s (n-l,E-e a )+Z s (n,E) 
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Here, the partition function is introduced, Z s (n, E) — 
^2 k F (Ek — E) , where the summation is taken over all states 
of n particles with the orbital s excluded. Correspondingly, 
the sum in Z a (n — 1, E — e s ) is taken over the states of n — 1 
particles with the orbital s excluded; this sum appears from 
the terms for which the orbital s is filled (n s = 1). For 
such states one can write Ek(n) — e a + Ek(n — 1) where 
Ek(n — 1) is the energy of the basis state with n — 1 par- 



ticles and e a — e a + 
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Note that to add the energy 



e a to Ek(n — 1) is the same as to subtract it from E because 
F = F(E k + e s - E). 

By taking e s independent on fc we assume the averaging 
over the basis states near the energy E is possible. Number 
of terms N in the partition function Z s is exponentially large, 
N = t — m! ,, , , therefore, one should consider In Z a which has 
slow dependence on n . In the case of large number of particles 
one can get In Z a (n— An, E — e a ) — In Z s (n, E) — a a An — (3 a l s 
where ot a = (3 a = ^g^S An = 1 . This leads to 

the "FD" type of the distribution, 



n a = (1 + exp(a a + /3 a i a ))~ 



(7) 



If the number of substantially occupied orbitals in the defini- 
tion of Z a is large, the parameters a a and /3 a are not sensitive 
as to which particular orbital s is excluded from the sum and 
one can assume a a — a = —u/T, (3 S — (3 = 1/T . The chem- 
ical potential {i and temperature T can be found from the 
conditions 
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In the case of many non- interacting particles (ideal gas) 
similar procedure transforms the canonical distribution (El) to 
the FD distribution (see e.g. ||). It is easy to check that 
the canonical distribution coincides with the FD distribution 
with a high accuracy even for very small number of particles 
provided the number of effectively occupied orbitals is large 
(when T > e or /j > e) . However, for the fixed total energy 
E, the temperature T in these two distributions is different. 
This difference can be explained by the dependence of the 
parameter a a on e s (see Eq. (Q)). Indeed, using expansion 
ot s ~ ol{ef) + ct'(e a — ef) one can obtain the relation between 
the FD (Pfd) and canonical (/3) inverse temperatures: (3fd = 
(3 + a'ep. 

Now, we can study the accuracy of the distributions (§,^,0) 
in the description of realistic quantum systems. As is known, 
the FD distribution is valid if the gas of particles can be 
considered as ideal, therefore, when the residual interaction 



is small enough, r <C u ~ ndo or T<C T (here do is the 
mean energy spacing for single-particle states). It is also as- 
sumed that the number of particles is large, n> 1. However, 
the equilibrium distribution for occupation numbers arises for 
much weaker condition, namely, when the number of princi- 
pal components N pc in exact eigenstates (see (0)) is large, 



f/D> 1 . In this case the fluctuations of the occu- 



pation numbers ~ N t 



-1/2 



Since the energy interval D for 
many-body states is exponentially small, it is enough to have 
relatively weak residual interaction Vb 3> D ~ do exp(— n) in 
order to get the equilibrium distribution. 

There are four regions of parameters depending on the 
strength of interaction and number of particles: 

(I) "regular" states, N pc ~ 1 for Vb < D ; 

(II) "initial chaotization" which is characterized by a rela- 
tively large number of "random" principal components, say, 
Np C ~ T/D > 10, however, the fluctuations are still large, 
iV- 1/2 > 0.1 ; 

(III) equilibrium F— distribution (Q) which is characterized 
by small fluctuations, N~ c 1/2 < 1 , or, N pc ~ F/D > 100; 
in this case components of eigenstates can be treated as ran- 
dom variables with the variance F and the Eq.(^j) gives actual 
distribution of the occupation numbers in quantum systems 
with interacting particles; 

(IV) canonical distribution (B which arises in the case of 
equilibrium plus large number of particles. If, in addition, the 
condition V <§C ndo is fulfilled, the standard FD distribution 
is valid. 

In practice, the condition (IV) of "thermalization" is not 
easy to satisfy in realistic systems like atoms or nuclei since 
n in this estimate is, in fact, the number of "active" particles 
(number of particles in the valence shell) rather than the total 
number of particles. Thus, the equilibrium F— distribution 
(||) which does not require the thermalization condition (IV) 
is more accurate. 

To test the above statements we have performed a detailed 
numerical study of the model of two-body random interac- 
tion. This model (see details in [B], [jjj) is described by few 
parameters: number n of particles, number m of orbitals and 
ratio Vb/do of the two-body interaction strength to the spac- 
ing between single-particle levels. For very small interaction 
the eigenstates are "regular" and the occupation number dis- 
tribution n s is a strongly fluctuating function even after av- 
eraging over a number of close eigenstates, see Fig.l. With 
an increase of the interaction keeping the number of parti- 
cles small, we obtain equilibrium F— distribution which is 
different from the FD distribution (see below). 

If, instead, we increase the number of particles keeping the 
interaction small, Vb <^ do , the distribution (^|) tends to the 
FD one (Fig. 2a). Finally, when the strength of the interac- 
tion is beyond the ideal gas approximation, the equilibrium 
distribution strongly deviates from the standard FD distri- 
bution (Fig. 2b). In fact, the latter result happens for a 
relatively small interaction Vb ~ O.ldo which, however, re- 
sults in the large value of spreading width, T > do since F 
increases with the number of particles very fast. 

However, the accuracy of the FD distribution can be im- 
proved by renormalization of the temperature. The point 
is that statistical effects of the interaction can be, at least 
in part, described by the increase of effective temperature. 
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For the first time, it was pointed out in Ref. |4j where we 
demonstrated that the spreading widths 7 of single-particle 
orbitals are imitated by the increase of temperature, namely, 
n s (e s , 7 ,T) « n s (e a ,j = 0,T + AT) = n s (e a ,f); for 7 < ^ 
one can obtain AT w 7 2 /16T. Below we present stronger 
arguments. Since we are mainly interested in the eigenstates 
which are in the lower half of the spectrum, the mean field 
energy (^) is substantially higher than the actual energy of 
eigenstates E^ . This results from level repulsion which forces 
lower energies to shift down. The mean field energy does not 
"know" about the level repulsion (the shift appears in the sec- 
ond order of perturbation theory), therefore, the unperturbed 
energy of the level with the same number would be a better 
estimate for it (E = Hu instead of E = E^'). Thus, we 
should substitute into Eq. (^|) the energy E which is higher 
than the energy of eigenstate: E = E^ + A^ ~ Hu. 

There is one more effect which leads to an additional shift 
A% . Since the density of states rapidly increases with the en- 
ergy, the number of admixed higher basis states is larger than 
that of the lower states (an extreme example is the ground 
state which contains the admixture of the higher basis compo- 
nents only). It is easy to estimate the corresponding increase 
of the mean field energy due to this effect if the spreading 
function F(Ek — E^ 1 ') is symmetric (to separate the effects 
we set here A^ — 0): 



A« = / F^(E k - E^)p (E k )dE h K d -^ffl(AEf (9) 



Here po is the unperturbed level density, and (AE) 2 is given 
by Eq. (^) . According to jjj , the shape of density of states for 
m > n > 1 is close to the Gaussian both for non-interacting 
and interacting particles, with the center of the spectrum E c 
and variances (ce) 2 and (5"b) 2 = (ce) 2 + (AE) 2 correspond- 
ingly. Therefore, one can obtain simple estimates of the shifts: 
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Taking into account these energy shifts and Eq.(g) one can 
estimate the increase of the effective temperature of the Fermi 
gas due to interaction. An accurate calculation of this effect 
requires more detailed knowledge of the spreading function F 
including it's weak "non-resonant" energy dependence on E k 
and E^\ However, the occupation numbers determined by 
the ^—distribution (^) are not sensitive to a particular choice 
of the spreading function provided the condition */ N pc S> I is 
fulfilled. To check this statement we have considered the form 
of the spreading function which takes into account important 
features of actual distribution F(E) . If the interaction is 
small, in the region not very far from the maximum, F max ~ 
D/F oc p _1 , the spreading function f(E) can be described 
by the Breit-Wigner form with the spreading width Few = 
2-kV 2 / D |J. Also, there is a shift of the maximum due to 
the repulsion between neighboring levels. These arguments 
allow us to suggest the improved expression for the spreading 
function valid for Yew < T: 
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where E = E {i) + A[ l) , Ti = Few and F 2 = F 2 /Ti. The 
value of F2 is found from the exact relation ([]) . Note that Fjf* 
from Eq. (0) automatically satisfies another exact relation 
y\ E^Fff' = Ek ; note, that the contribution of the energy 
shift A^ is compensated by the increase of the level density 

p (eW). ' 

It follows from Eq.( lllh that the shift of energy due to the 
enhanced admixture of higher basis components is two times 
smaller than in Eq. ^ since the increase of density po is 
partly compensated by the factor po~ x ^ 2 in Eq.([lf|): 



dE 



(AE) 2 ~ A^ 
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Thus we should substitute to Eq.(||) the value E = E (i) + 
A^ + A { 2 i} ~ 2Hu - E (i) . After taking into account 
this shift, the FD expression gives the same result as the 
F— distribution; namely, in Fig. 2b the FD curve (circles) is 
shifted to that given by stars, both curves also agree with 
numerical experiment (Fig. 3). 

For small Vb and large number of particles we have A[ <§C 
ri <g; F2 , therefore, in the central part the distribution ( |ll| ) 
has the Breit-Wigner shape with the width Few ■ However, 
numerical calculations 0] H demonstrate that at larger 
interaction Vb the width of the spreading function becomes 
linear in Vb and it is better to put Fi = F2 = F. One can 
write the extrapolation expression both for small and large 
values of V (see also ||) Ti = F B wF/(F B w + F). We have 
also checked that for large number N pc of principal compo- 
nents the distribution of the occupation numbers does not 
depend on Ti; this fact can be treated as a signature of the 
equilibrium. 

Finally, we discuss the transition to mesoscopic systems. 
The result depends on the dimensionality d of the system since 
d ~ l~ 2 ; Vb ~ l~ d , therefore, Vo/do ~ r (d_2) where I is the 
size of the system. Thus, for d — 1 one has Vb S> do which 
means that strong mixing (chaos) starts just from the ground 
state. This is in accordance with the absence of a gap in 
the distribution of occupation numbers in ID case (Luttinger 
liquid). In the 3D case we have Vb <C do and the admixture 
of the higher states to the ground state can be considered 
perturbatively which is consistent with the non-zero gap at 
T — . One can see that the transition between regular region 
(I) and the initial chaos region II in the 3D case occurs for 
high states when Vb > D ~ do exp(— n|) where n" 2 oc E 1 ^ 2 is a 
number of excited particles. On the other hand, the transition 
from the fluctuating to the equilibrium regime, (II)— > (III), 
requires the decrease of D only by one order of magnitude: 
n 3 ~ n 2 + 2 . This means that such transition in 3D is a very 
sharp "phase transition". 

In conclusion, we developed new method for the calcula- 
tion of the occupation numbers in finite systems of interacting 
particles. The method is based on the assumption that ex- 
act eigenstates are "chaotic" superpositions of the shell-model 
basis states, and the smooth spreading function for the eigen- 
states components can be introduced. This assumption and 
the results are confirmed by numerical experiments. The "mi- 
crocanonical" partition function which we have introduces can 
be used for further studies of statistical and thermodynamic 
properties of finite systems of interacting particles. 
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We also demonstrated that occupation numbers in the sys- 
tems of interacting particles can be reasonably described by 
Fermi-Dirac distribution with renormalized parameters. As 
usual, mean field (and possibly other "regular" effects) can 
be included into single-particle energies e s . Staistical effects 
of the residual interaction (mainly due to non-diagonal ma- 
trix elements of the Hamiltonian matrix) increase effective 
temperature. 
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for the hospitality during his visit when this work was done. 
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(b) strong interaction, Vo = 0.08, F w 16.6 , N pc ^ 11000. 

Fig. 3 Distribution of occupation numbers for Ce atom pa- 
rameters n = 4, m = 11, do = 1,Vq = 0.12, N pc ~ 48 (boxes) 
in comparison with the FD distribution, with the corrected 
temperature (circles) and with the F— distribution (stars). 
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Fig.l. Distribution of occupation numbers n s for n = 4 
particles and m = 11 orbitals for the two-body random inter- 
action model 0] with do = 1 , weak interaction Vo = 0.04 
and total energy E = 17.33 (to compare with the ground 
state energy Ef ~ 12.1) . Dashed boxes represent numeri- 
cal data averaged over 20 Hamiltonian matrices of the size 
N = 330 with different realization of the random interaction 
and over the energy range E ± SE with SE — 0.25. Circles 
correspond to the FD distribution with temperature T and 
chemical potential /i defined by total number n of particles 
and total energy E . Stars are given by the F— distribution 
(^) for n s (E) . The number of principal components in exact 
eigenstates is N pc ~ 13.0 . The latter has been calculated via 
the inverse participation ratio, N pc = (^2 k (Ek ~ E)) -1 . 
Large fluctuations ~ N pc are seen. 

Fig. 2. Distribution of occupation numbers n s for larger 
number n = 14 of particles and m = 28 of orbitals with 
the same do and temperature. Direct diagonalization of huge 
Hamiltonian matrices is not possible in this case. Circles are 
the FD distribution; stars are the F— distribution (|^). 

(a) weak interaction, V = 0.003, F « 0.62 , N pc » 60. 
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